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The motion of molecules on solid surfaces is of interest for technological applications such as 
catalysis and lubrication, but it is also a theoretical challenge at a more fundamental level. The 
concept of activation barriers is very convenient for the interpretation of experiments and as input 
for Monte Carlo simulations but may become inadequate when mismatch with the substrate and 
molecular vibrations are considered. We study the simplest objects diffusing on a substrate at 
finite temperature T, namely an adatom and a diatomic molecule (dimer), using the Langevin 
approach. In the driven case, we analyse the characteristic curves, comparing the motion for different 
values of the intramolecular spacing, both for T — and T ^ 0. The mobility of the dimer 
is higher than that of the monomer when the drift velocity is less than the natural stretching 
frequency. The role of intramolecular excitations is crucial in this respect. In the undriven case, the 
diffusive dynamics is considered as a function of temperature. Contrary to atomic diffusion, for the 
dimer it is not possible to define a single, temperature independent, activation barrier. Our results 
suggest that vibrations can account for drastic variations of the activation barrier. This reveals a 
complex behaviour determined by the interplay between vibrations and a temperature dependent 
intramolecular equilibrium length. 



I. INTRODUCTION 

The diffusion of adatoms and molecules on surfaces is 
recently attracting much attention, in order to under- 
stand many properties of technological interest, such as 
thin film growth, catalysis and dissociation [1,2]. Useful 
experimental techniques for probing the surface diffusion 
of molecular adsorbates have been developed. These al- 
low to follow the dynamical details of surface diffusion, 
using Scanning Tunneling Microscopy [3], and to im- 
age physisorbed atoms and clusters, using Field Ion Mi- 
croscopy [4] . Most theoretical works have focussed on the 
determination of energy barriers for diffusion in different 
systems [5-9] , usually on the basis of energy arguments 
and neglecting the role of internal degrees of freedom. 
However, it is often suggested that diffusion dynamics 
can be strongly affected by the presence of intramolec- 
ular motion [10,11]. Here we address this problem for 
the simple but important case of a dimer physisorbed 
on a periodic substrate, trying to link the macroscopic 
diffusive behaviour to microscopic degrees of freedom. 
Our theoretical model contains the salient features of a 
diatomic molecule physisorbed on a periodic substrate. 
However, our work does not aim to probe diffusion in a 
specific molecule-surface system, but to understand dif- 
fusion mechanisms expected for systems of this kind. 

In Sec. II we give some details of the model we have 
used. Sec. Ill is devoted to the discussion of the results in 
the presence of an external driving, while Sec. IV deals 
with pure thermal diffusion. Some concluding remarks 
are presented in the last section. 



II. MODEL 

We consider a monomer (adatom) and a dimer (di- 
atomic molecule) moving on a periodic one-dimensional 
substrate at finite temperature. The particle-substrate 
interaction is modelled by a periodic function U. Specif- 
ically, for the monomer 

U{x) = U (l - cos(kx)) (1) 

and for the dimer 

U(xi, X2) — Uq(2 — cos(fca;i) — cos(fca;2)), (2) 

where x represents the spatial coordinate and fc = 2ir/a, 
a being the lattice constant of the substrate. 

For the dimer the interparticle interaction is taken to 
be harmonic: 

K 

V(x!,x 2 ) = yfe ~ %i ~ I) 2 , (3) 

where / is the spring equilibrium length. In order to take 
into account the finite temperature T of the substrate, we 
adopt the Langevin approach. In this framework the mo- 
tion of the sliding object is described exactly, while the 
substrate is treated as a thermal bath. The equations of 
motion of the monomer and the dimer are respectively 

mx + mrjx = —Uq sin(fcx) + f + F (4) 
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and 



mx\ + mr)Xi 
mx 2 + mr]x 2 



K(x 2 - Xl -l)-U sm(k Xl ) + .h + F 
K(x! -x 2 + l)-U sm{kx 2 ) + h + F, 

(5) 



where an explicit damping term mr/ii modelling energy 
dissipation has been introduced (77 is the phenomenolog- 
ical friction coefficient), and the effect of finite temper- 
ature T is taken into account by the stochastically fluc- 
tuating forces /j. These two terms are related by the 
dissipation-fluctuation theorem: 



</i(*)/j(0) >=2mk B T5ij5(t). 



(6) 



external force F. For the monomer, we can rely on the re- 
sults of Kramers' Transition State Theory (TST) [12-14], 
yielding analytical results for certain regimes of friction 
and for relatively high values of the diffusion barrier 2Uq. 
We present in Fig. 1 the results obtained with our simu- 
lation for Uq = 2.5 and r\ = 1, which we will use as refer- 
ence in discussing the behaviour of the dimer. Fig. 1(a) 
shows the velocity-force characteristic in the T = case. 
Increasing F adiabatically the velocity remains zero un- 
til a critical force F™ on is reached, since the system has 
to overcome a finite static friction force in order to get 
out of the potential well. Then the particle slides with 
constant velocity v on the substrate, resulting in a linear 
relation between v and F for high values of the driving 
force: 



It is convenient to rewrite the equation of motion in adi- 
mensional units by introducing a characteristic time 



F ~ i]v. 



(10) 



\k B Tk 2 



1/2 



and defining 

x = kx, t = t/r, fj = rjT, Uq = Uo/iksT) 



f = f/(kk B T), F = F/(kk B T), l = kl, 



K = K/(k 2 k B T). 3 



For typical values m ~ 2 • lCT 26 ^, T = 300K, a ~ 2l, 
we have r ~ 0.25ps. In this way Eqs. (4) and (5) become 
(in the following we omit the tildes for simplicity) 



x + r\x = — t/o sin a; + / + F 



(7) 
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and 




(xi+r/xx = K(x 2 - xi - I) - U sin xi + fi + F 
\ x 2 + r)± 2 = K(xi - x 2 + I) - Uq sin x 2 + f 2 + F, 

(8) 

and the fluctuation-dissipation relation Eq. (6) 

<fi(t)f j (0)>=2riSi j 5{t). (9) 

We perform Molecular Dynamics (MD) simulations, in- 
tegrating the equations of motion using a velocity- Verlet 
algorithm, with time step A = 1CP 4 t, averaging the tra- 
jectories over several hundreds of realizations (~ 300 in 
driven case and ~ 3000 in the undriven case with F = 0) , 
in order to reduce the statistical noise due to the stochas- 
tic term. 

III. DRIVEN CASE 

We have analyzed the response properties of Eqs. (7) 
and (8) by studying the characteristic curves, i.e. the 
behaviour of the mobility of the system subjected to the 
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FIG. 1. Characteristic curve of the monomer for T = (a) 
and r / (b). The arrows in (a) indicate the parts of the 
curve where the force is increased or decreased. The lines with 
crosses in (a) and (b) are the numerical data, while the dashed 
line and the dotted line in (b) are fits to the data respectively 
according to Eq. (11) and Eq. (12). The parameters used are 
Uo = 2.5 and r] = 1. 

If we decrease the force adiabatically a hysteresis effect 
is observed, since the mobility is different from zero also 
when F < F^ on and vanishes at F = F 2 mo ™ < F^ on . 
This corresponds to a bistability between the locked and 
the running solution, which is present when the friction 
coefficient i] is not very large, t]/\/Uq < 1.19 [14]. In 
the case T ^ this bistability disappears, at least for 
intermediate-large values of r/. The system has always 
the chance to get untrapped because of thermal fluctua- 
tions. As a consequence the mobility is always different 
from zero when F ^ 0, as shown in Fig. 1(b). If 77 is 
relatively large, it is possible to describe analytically the 
behaviour for small F by means of Kramers' TST [12]. 
The result for the behaviour of v(F) is 
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+ Uq cos arcsin 



1/2 



F 

exp(-2C/ ) (exp(Fa/2) - exp(-Fa/2)) (11) 



As we can see from Fig. 1(b), Eq. (11) reproduces the 
simulation results for F < 1, but this range increases 
for larger values of Uo, where the activated processes 
are more pronounced. The behaviour of the mobility for 
large F, where the particle performs a drift motion with 
a small contribution of the noise term can be shown to 
be [15] 



(12) 



In particular, if F is very large compared to [To the term 
in parenthesis tends to one and the characteristic curve 
becomes linear: 



F ~ rjv. 



(13) 



Now we consider the driven dimcr described by Eq. (8). 
In this case, the presence of the interparticle interac- 
tion renders the problem more complex, giving rise to a 
richer dynamical behaviour. We have analyzed the char- 
acteristic curves for different values of the intramolecular 



spacing /, namely I 



I = a/2 and I — r g a where 



T g = (1 + V5)/2 is the golden mean, both for T = and 
T^O (see Fig. 2), and compared them to the monomer 
characteristics (see Fig. 3). 



For the dimer it is convenient to rewrite the equation 
of motion using the centre of mass (CM) and relative 
coordinates, defined by 



xcm = (xi + x 2 )/2 
From Eq. (8) we obtain 



Vr 



X2 — X\ — 1 



(14) 



'icM = -tjxcm - U cos((y r + l)/2) smxcM + F 
jjr = -VVr - 2Ky r - 2U sm((y r + l)/2) cosx C m 



(15) 



In our MD simulations we choose the initial configura- 
tion which minimizes the total potential energy. As for 
the monomer a critical force Ff lm , which depends on the 
value of /, is needed to achieve motion for T = 0. Then 
for larger values of F the velocity increases as a func- 
tion of the external force, but at a certain value of the 
force F$ lm another plateau appears in the vcm — F plane 
( v c m being the CM mean velocity) , signalling a dynami- 
cal crossover in the system. Finally, keeping on increasing 
the force, the linear regime is recovered (Fig. 2(a)). In 
the CM frame, the external potential leads, for a drift 
motion xqm ~ vcutt, to a time-periodic force acting on 
the particles, with "washboard" frequency given by vcm- 
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FIG. 2. Characteristic curves of the dimer for T — 
(a),(b) and T / ( c ),(d) for three values of I. The vertical 
dot-dashed line passes through vcm = <^o- The parameters 
used are Uo = 2.5, T) = 1 and K = 1. 

Note that we have chosen I — r g a because it is a pro- 
totypical example of an incommensurate length ratio. 
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FIG. 3. Comparison between the characteristic curves of 
the monomer and the dimer for T = (a) and T / (b). 
The parameters used are the same as in Figs. 1 and 2. 

The force F^ tm , where the second plateau appears, 
physically corresponds to the point where the washboard 
frequency is in resonance with the stretching frequency 
of the dimer ix>o = {2K) 1 / 2 , exciting the internal degrees 
of freedom. This resonance mechanism was also found in 
the Frcnkel-Kontorova model in the low friction limit [16] . 
For the dimer two regions of bistability can be observed, 
corresponding to the two plateaus in the characteristic 
curve; thus, when F is decreased a first hysteresis oc- 
curs in proximity of vcm = ^0, giving a critical value 
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pdim pdim ^ w h ere ^he characteristic curve has a dis- 
continuous derivative, while a second hysteresis is found 
when F is decreased further from Ff lm and another crit- 
ical value F% lm < Ff lm , where v C m = 0, is obtained 
as for the monomer (these hysteresis curves are shown 
only for the / = r g a dimer in Fig. 2). For a closer com- 
parison between the monomer and the dimer at T = 
see Fig. 3(a). We note that the qualitative behaviour is 
the same for different I but the values of F f tm and F^" 1 
can differ significantly as a function of I. This is due 
to energetic reasons: the dimer with I = a/2 is energet- 
ically favourite since, on average, it has to overcome a 
lower barrier, while for I = a the two particles tend to 
be pinned in the minima and to behave like a monomer 
(in fact the characteristic curve of the I = a dimer is 
practically superimposed on that of the monomer). The 
resonance frequency is highlighted in Fig. 2(b) by plot- 
ting F — 7]vcm as a function of vcm- a clear peak of the 
characteristic curves at vcm — <^o is visible when I ^ a. 



In order to understand better the resonance mecha- 
nisms we consider the case / — a at T = 0. Assuming 
a CM drift motion, xcm(^) = + ucm*, and lineariz- 
ing the equation of motion (15) for I = a we obtain the 
following equation for y r : 
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FIG. 4. Relative and CM motion of the dimer for / = a 
with initial condition xi — x\ = 1.27a, for different values of 
the external driving F. y r is shown in (a) and (c), while the 
deviation of xcm — xo from vcmt is plotted in (b) and (d) 
(xo is the initial condition for CM). The parameters used are 
the same as in Fig. 2. 
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+ r\y r + 2Ky r = U cos(v C Mt + x )y r . (16) 



This is the equation of a parametric oscillator, for which 
an exponential increase of the amplitude is expected for 
vcm — 2wo within a given instability window, which we 
estimate to be 2.21 < vcm < 2.37. We note that indeed 
the amplitude of y r increases exponentially in a certain 
range of F, as shown in Fig. 4, but it saturates at long 
times. This is due to the fact that Eq. (16) assumes a con- 
stant CM velocity. Actually, in the full system Eq. (15) 
xcm is coupled to y r , so that vcm decreases slightly dur- 
ing the dynamics as shown in Fig. 4(d). This is enough 
to shift vcm out of the instability window, thus stopping 
the increase of y r . 



The behaviour outlined for T = smears out at fi- 
nite temperatures. Although a hysteretic behaviour has 
been found for long periodic chains in the low friction 
limit [17], no bistabilities and hysteresis are present in our 
shown in Fig. 2(c). The static friction force van- 
ishes and the mobility is sensitive to the value of I. In par- 
ticular the dimer with I = a has the lowest mobility. Note 
that the curves for different I still cross at vcm = ^o- 
This crossover of the mobility can be observed also com- 
paring the monomer and the dimer (Fig. 3(b)). More- 
over, as it can be seen from Fig. 2(d), plotting F — tjvcm 
vs. vcm, the resonance peak at vcm — still survives 
for I 7^ a. For I = a the peak is shifted to a smaller value 
of vcm- This is due to the fact that in the commen- 
surate case, even though the static friction force found 
at T — vanishes, a larger force is needed to reach the 
sliding state and this corresponds to the point of highest 
curvature in the characteristic plotted in Fig. 2(c). To a 
certain extent this resembles the monomer case, where a 
very similar behaviour is found. 



This shows how intramolecular vibrations can be reso- 
nantly excited due to the sliding on a periodic substrate 
and that the details of the resulting relative motion arc 
non trivial. For instance, whether this could represent 
a mechanism for dissociation depends on the maximum 
excursion from the equilibrium distance. Moreover, the 
very nature of parametric resonances makes the tempo- 
ral behaviour very much dependent on the initial values 
of the interatomic spacing which is in turn related to vi- 
brational energy and temperature. 



IV. DIFFUSION 

Next we examine the pure thermal diffusion, i.e. we 
study the motion of the monomer and the dimer with 
F = 0. This problem is computationally more time con- 
suming, since the motion is completely random and an 
averaging over many MD realizations is needed. We ex- 
tract information about the diffusive behaviour of the 
particles by computing the Mean Square Displacement 
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(MSD) < x 2 (t) >, where < • > denotes the average over 
the realizations. The diffusion coefficient D is defined by 



D 



lim 

t— *oo 



< X 2 {t) > 

2t ' 



(17) 



Usually it is assumed that the dependence of D on tem- 
perature should follow the Arrhenius law: 



D = D exp{-E a ) 



(18) 



where Do is a prefactor and E a is the activation energy 
for diffusion, scaled to fc^T. In Eq. (18) both D and E a 
are assumed to be T independent. However, some recent 
studies have already shown that there may be deviations 
from the Arrhenius behaviour [10,11,18]. Fig. 5 shows 
the diffusion coefficient D as a function of the single- 
particle energy barrier Eb = 2U~o (scaled to fc^T) for the 
monomer and the dimer for different values of I and two 
values of K. We note that while the monomer has a 
uniquely defined activation energy E a ~ 0.95, a clear de- 
viation from the Arrhenius law is observed in the dimer 
case, except for the case I = a, where it is less evident (as 
explained in Sec. Ill the commensurate dimer dynamics 
is similar to that of the monomer). 



10 



0.1 



0.01 







\\ \ 




V\ 




k \ 




Hi v * 




monomer 
dimer l=a/2 (K=1) 
dimer l=a (K=1) 
dimer l=xa (K=1 ) 
dimer l=xa (K=0.1 ) 


---*< — 
— sie- - - - 

El 



FIG. 5. Diffusion coefficient of the monomer and the dimer 
for different values of I and K vs. Et = 2Uo- Here F = 
and T] — 0.1. Note that for the free diffusion case Eb —> the 

diffusion coefficient of the monomer is - = 10 while that of 

i 

the dimer is ^ — 5 (independent of I). 



(low Ei,). A similar mechanism was also proposed in a 
model of heteroepitaxial island diffusion [11]. Since we 
are considering small values of E\, (rigorously Eq. (18) 
should be valid for Eb — > oo), finite-barrier effects are also 
important, as illustrated in Ref. [18]. We are currently 
trying to unravel the latter from the effects due to vibra- 
tions. Our results suggest that the role of intramolecular 
vibrations is relevant for the diffusion dynamics. The 
possibility to excite internal degrees of freedom, even for 
a simple diatomic molecule, can affect the temperature 
dependence of the activation energy [10] and determine 
a complex diffusion behaviour which can depend on the 
lattice commensurability and on the interatomic spacing. 
We plan to address this issue more extensively in a future 
work. 



V. CONCLUSIONS AND PERSPECTIVES 

We have presented a simplified model to study the 
driven and undriven motion of monomers and dimers 
on a periodic substrate. We have pointed out the pe- 
culiar effects of the substrate on the CM and relative 
motion of the dimer, such as nonlinear mobility, bistabil- 
ities and resonance processes. In this respect, the cou- 
pling between the internal degrees of freedom is crucial. 
In particular, the diffusion dynamics reveals strikingly 
complex features determined both by energetic mecha- 
nisms and by the role of intramolecular motion. It would 
be worthwhile to explore further this issue going from 
the simple ID model to a more realistic approach, con- 
sidering for example the motion and the orientation of 
large molecules on a 2D surface and the diffusion of long 
chains. 
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The I — a/2 dimer diffuses faster than the monomer 
(at least for Eb > 1) due to the most favourable energy 
configuration (sec also [6]). Note that the crossover point 
in the diffusion behaviour can depend on the force con- 
stant K, as it can be seen by comparing the two curves 
for I = r g a with K = 1 and K = 0.1. The dependence of 
the activation energy on temperature (or equivalently on 
the energy barrier) can be attributed to a temperature 
dependent intramolecular length, which affects the diffu- 
sion behaviour especially in the high temperature regime 
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